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ABSTRACT 



Context. Many chemical changes occur during the collapse of a molecular cloud to form a low-mass star and the surrounding disk. 
One-dimensional models have been used so far to analyse these chemical processes, but they cannot properly describe the incorpora- 
tion of material into disks. 

Aims. The goal of this work is to understand how material changes chemically as it is transported from the cloud to the star and the 
disk. Of special interest is the chemical history of the material in the disk at the end of the collapse. 

Methods. A two-dimensional, semi-analytical model is presented that follows, for the first time, the chemical evolution from the pre- 
stellar core to the protostar and circumstellar disk. The model computes infall trajectories from any point in the cloud and tracks the 
radial and vertical motion of material in the viscously evolving disk. It includes a full time-dependent radiative transfer treatment of 
the dust temperature, which controls much of the chemistry. A small parameter grid is explored to understand the effects of the sound 
speed and the mass and rotation of the cloud. The freeze-out and evaporation of carbon monoxide (CO) and water (H2O), as well as 
the potential for forming complex organic molecules in ices, are considered as important first steps to illustrate the full chemistry. 
Results. Both species freeze out towards the centre before the collapse begins. Pure CO ice evaporates during the infall phase and 
re-adsorbs in those parts of the disk that cool below the CO desorption temperature of ~18 K. H2O remains solid almost everywhere 
during the infall and disk formation phases and evaporates within ~10 AU of the star. Mixed CO-H2O ices are important in keeping 
some solid CO above 18 K and in explaining the presence of CO in comets. Material that ends up in the planet- and comet-forming 
zones of the disk (~5-30 AU from the star) is predicted to spend enough time in a warm zone (several lO'' yr at a dust temperature 
of 20--40 K) during the collapse to form first-generation complex organic species on the grains. The dynamical timescales in the hot 
inner envelope (hot core or hot corino) are too short for abundant formation of second-generation molecules by high-temperature 
gas-phase chemistry. 

Key words, astrochemistry - stars: formation - circumstellar matter - planetary systems: protoplanetary disks - molecular processes 



1. Introduction 

The formation of low-mass stars and their planetary systems is 
a complex event, spanning several orders of magnitude in tem- 
poral and spatial scales, and involving a wide variety of physical 
^ and chemical processes. Tha nks to observ ations (see reviews by 
' Idi Francesco et alJlIOOTi and lWhite et al.1 12007). theory (see re- 
vi ew bvlshu et al.lll987h and computer sim ulations (see reviews 
bv lKlein et alj2007l and lDullemorid et al.l2 007). the general pic- 
ture of low-mass star formation is now understood. An insta- 
bility in a cold molecular cloud leads to gravitational collapse. 
Rotation and magnetic fields cause a flattened density structure 
early on, which evolves into a circumstellar disk at later times. 
The protostar continues to accrete matter from the disk and the 
remnant envelope, while also expelling matter in a bipolar pat- 
tern. Grain growth in the disk eventually leads to the formation 
of planets, and as the remaining dust and gas disappear, a mature 
solar system emerges. While there has been ample discussion in 
the literatur e on the or igin and evolution of grains in disks (see 
reviews by iNatta et al]i200'7i and iDominik et alj|2007l) . little at- 
tention has so far been paid to the chemical history of the more 
volatile material in a two- or three-dimensional setting. 

Chemical models are required to understand the observations 
and develop the simulations (see reviews by ICeccarelli et alj 



l2007llBergin et alJl2007l and lBergin & Tafallall2007h . The chem- 
istry in pre-stellar cores is relatively easy to model, because 
the dynamics and the temperature structure are simpler before 
the protostar is formed than afterwards. A key result from the 
pre-stellar core models is the depletion of many carbon-bearing 
species toward s the centre of the core dBergin & Langeij|l997l : 
iLee et alJl2004l). 

ICeccarelli et al.l (Il996h modelled the chemistry in the 
collapse phase, and others have done so more recently 
(Rodgers & Charnlev 2003; Doty et al. 2004; Lee et al. 2004, 
Garrod & Herbsti,2006; .Aikawa et al..,2008i; ,Garrod et al.il2008l) . 
All of these models are one-dimensional, and thus necessarily 
ignore the circumstellar disk. As the protostar turns on and heats 
up the surrounding material, all models agree that frozen-out 
species return to the gas phase if the dust temperature surpasses 
their evaporation temperature. The higher temperatures can fur- 
ther drive a hot-core-like chemistry, and complex molecules may 
be formed if the infall timescales are long enough. 

If the model is expanded into a second dimension and the 
disk is included, the system gains a large reservoir where in- 
falling material from the cloud can be stored for a long time 
before accreting onto the star This can lead to further chem- 
ical enrichment, especially in the warmer parts of the disk 
jAikawa et al.ll 1 9971; IXikawa & Herbstll 1 9991; IWiUacv & Langerl 
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Aikawa et al.i i2008 l). The interior of the disk is shielded from 



direct irradiation by the star, so it is colder than the disk's sur- 
face and the remnant cloud. Hence, molecules that evaporated as 
they fell in towards the star may freeze out again when they en- 
ter the disk. This was first shown quantitatively bv lBrinch et al.l 
( I2OO 8', hereafter BWH08) using a two-dimensional hydrody- 
namical simulation. 

In addition to observations of nearby star-forming regions, 
the comets in our own solar system provide a unique probe 
into the chemistry that takes place during star and planet for- 
mation. The bulk composition of the cometary nuclei is be- 
lieved to be largely pri stine, closely reflecting the c omposition 
of the pre-solar nebula (IBockelee-Morvan et al.l2004b . However, 
large abundance variations have been observed between i ndivid- 
ual c omets and these remain poorly understood (Kobayas hi et all 
12007 ). Two-dimensional chemical models may shed light on the 
cometary chemical diversity. 

Two molecules of great astrophysical interest are carbon 
monoxide (CO) and water (H2O). They are the main reservoirs 
of carbon and oxygen and control much of the chemistry. CO is 
an important precursor for more complex molecules; for exam- 
ple, solid CO can be hydrogenated to forma ldehyde (H2CO) and 
methanol (CH3OH) at low temperatures (IWatanabe & Kouchil 
12002 : Fuchs et al. subm.). In turn, these two molecules form 
the basis of even larger organic species like methyl fo rmate 
(HCOOCH3; lGarrod&Herbsil2006t iGarrod et al.ll2008h . The 
key role of H2O in the formation of life on Earth and potentially 
elsewhere is evident. If the entire formation process of low-mass 
stars and their planets is to be understood, a thorough under- 
standing of these two molecules is essential. 

This paper is the first in a series aiming to model the chem- 
ical evolution from the pre-stellar core to the disk phase in two 
dimensions, using a simplified, semi-analytical approach for the 
dynamics of the collapsing envelope and the disk, but includ- 
ing detailed radiative transfer for the temperature structure. The 
model follows individual parcels of material as they fall in from 
the cloud into the disk. The gaseous and solid abundances of CO 
and H2O are calculated for each infalling parcel to obtain global 
gas-ice profiles. The semi-analytical nature of the model allows 
for an easy exploration of physical parameters like the cloud's 
mass and rotation rate, or the effective sound speed. Tracing the 
temperature history of the infalling material provides a first clue 
into the formation of more complex species. The model also pro- 
vides some insight into the origin of the chemical diversity in 
comets. 

Section|2]contains a full description of the model. Results are 
presented in Section |3] and discussed in a broader astrophysical 
context in Section|4] Conclusions are drawn in Section|5] 

2. Model 



The physical part of our two-dimensional axisymmetric model 
describes the collapse of an initially spherical, isothermal, 
slowly rotating cloud to form a star and circumstella r disk. 
The collapse dynamics are taken from IShul d 19771 here- 
after S77 ), including the effects of rotation as described by 
ICassen & MoosmanI (Il98li hereafter CM81) and lTerebev et al.1 
(11984 hereafter TSC84). Infalling material hits the equatorial 
plane inside the centrifugal radius to form a disk, whose fur- 
ther evolution is constrained by conservation of angular mo- 
mentum (L ynden- Bell & Pringle 1974). So me prop e rties o f the 
star and the disk are adapted from 1 Adams & Shul ( Il986l) and 
lYoung & EvansI (l2005l hereafter YE05). Magnetic fields are not 



included in our model. They are unlikely to aff'ect the chemistry 
directly and their main physical effect (causing a flattened den- 
sity distribution; Galli & Shu 1993) is already accounted for by 
the rotation of the cloud. 

Our model is an extension of the one used by 

iDullemond et al.l (l2006h to study the crystallinity of dust 
in circumstellar disks. That model was purely one-dimensional; 
our model treats the disk more realistically as a two-dimensional 
structure. 



2.1. Envelope 

The cloud (or envelope) is taken to be a uniformly rotating singu- 
lar isothermal sphere at the onset of collapse. It has a solid-body 
rotation rate Q^) and an r"^ density profile (|S77|) : 



Poir) 



(1) 



where G is the gravitational constant and Cs the effective sound 
speed. Throughout this work, r is used for the spherical radius 
and R for the cylindrical radius. Setting the outer radius at renv, 
the total mass of the cloud is 



Mo = 



2ch 



(2) 



After the collapse is triggered at the centre, an expansion 
wave (or coll apse front) travels outwards at the sound speed 
(IsTTL ITSC841) . Material inside the expansion wave falls in to- 
wards the centre to form a protostar The infalling material is 
deflected towards the gravitational midplane by the cloud's rota- 
tion. It first hits the midplane inside the centrifu gal radius (where 
gravity balances angular momentum; ICM8II) . resulting in the 
formation of a circumstellar disk (Section lZ2b . 

The dynamics of a collapsing singular isothermal sphere 
were computed by S77 in terms of the non-dimensional variable 
X - r/cj, with t the time after the onset of collapse. In this self- 
similar description, the head of the expansion wave is always 
at X - I. The density and radial velocity are give n by the non- 
dimensional variables and v, respectively. (|S77| uses a for the 
density, but our model already uses that symbol for the viscosity 
in Section l272l ) These variables are dimensionalised through 



p(r, t) = 



AnOfi ' 



Ur{r, t) — CsV(x) . 



(3) 



(4) 



Val ues for and v are tabulated in lS7" 

|CM81 and TSC84 analysed the effects of slow uniform rota- 
tion on the S77 collapse solution, with the former focussing on 
the flow onto the protostar and the disk and the latter on what 
happens further out in the envelope. In the axisymmetric ItSCS^ 
description, the density and infall velocities depend on the time, 
f, the radius, r, and the polar angle, 8: 



Jl{x, 6, T) 
AnGfi 



Urir, 6, t) — Csv(x, 6, t) , 



(5) 



(6) 



where r = Qof is the non-dimensional time. The polar velocity 
is given by 



ueir, 6, t) — Csw{x, 6, t) . 



(7) 
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The differential equations from lTSC84l were solved numerically 
to obtain solutions for v and w. 

Th e [TSC84 solution breaks down around x - t^, so the 
ICMSlI solution is used inside of this point. A streamline through 
a point (r, 6) effectively originated at an angle Gq in this descrip- 
tion: 



cos 60 - cos 6 Rr 
— ^ =0, 



(8) 



sin" &o cos fo 

where is the centrifugal radius, 
1 

Ye 

with niQ a numerical factor equal to 0.975. The ICMSTl radial and 
polar velocity are 



(9) 



cos 6 



r V cos 00 



ue(r, 0, f) 



1 + 



IGM 
r 

and the ICM8II density is 
M 



cos 6 cos 00 - cos 
cos 00 sin 



p(r,0,f) = - 



Anr^Ur 



R t'^ 
1 +2— P2(cos0o) 
r 



(10) 



(11) 



(12) 



where P2 is the second-order Legendre polynomial and M - 
m()cl/G i s the total acc retion rate from the envelope onto the star 
and disk (IS77l:lTSC84l) . The primary accretion phase ends when 
the outer shell of the envelope reaches the star and disk. This 
point in time (t^^c - Mq/M) is essentially the beginning of the T 
Tauri or Herbig Ae/Be phase, but it does not yet correspond to a 
typical T Tauri o r H erbig A e/Be object (see Section lTSb . 

The lTSC84l and lCM8l1 solutions do not reproduce the cavi- 
ties created by the star's bipolar outflow, so they have to be put 
in separately. O utflows have been o bserved in two shapes: con- 
ical and curved dPadgett et al.lll999l) . Both can be characterized 
by the out flow op ening angle, 7, which grows with the age of 
the object. lArce & Sargent (2006) found a linear relationship in 
log-log space between the age of a sample of 17 young stellar ob- 
jects and their outflow opening angles. Some explanations exist 
for the outflow widening in general, but it is not yet understood 
how yit) depends on parameters like the initial cloud mass and 
the sound speed. It is likely that the angle depends on the relative 
age of the object rather than on the absolute age. 

The purpose of our model is not to include a detailed de- 
scription of the outflow cavity. Instead, the outflow is primar- 
ily included because of its effect on the temperature profiles 



dWhitnev et al 



2003h Its opening angle is based on the fit by 



lArce & Sargent l2006h to their Fig. 5, but it is taken to depend 



on f/facc rather than t alone. The outflow i s also kept smaller, 
which brings it closer to the IWhitnev et al.l angles. Its shape is 
taken to be conical. With the resulting formula. 



y(t) t 
log — = 1.5 + 0.26 log — . 

deg face 



(13) 



the opening angle is always 32° at f = face- The numbers in 
Eq. ( fTST l are poorly constrained; however, the details of the out- 
flow (both size and shape) do not affect the temperature pro- 
files strongly, so this introduces no major errors in the chem- 
istry results. The outflow cones are filled with a constant mass 
of O.OO2M0 at a uniform density, which decreases to 10"'-10'* 
cm~^ at 4cc depending on the model parameters. The outflow 
effectively removes about 1% of the initial envelope mass. 



2.2. Disk 

The rotation of the envelope causes the infalling material to be 
deflected towards the midplane, where it forms a circumstel- 
lar disk. The disk initially forms inside the centrifugal radius 
dCMBll) . but conservation of angular momentum quickly causes 
the disk to spread beyond this point. The evolution of the disk is 
governed by viscosity, for which our m odel uses the common a 
prescription ( Shakura & Sunvaevll973h . This gives the viscosity 
coefficient v as 



v(R, t) = ac^^H . 



(14) 



The sound speed in the disk, Cs,d = ^jkTml^m-p (with k the 
Boltzmann constant, njp the proton mass and yu the mean molec- 
ular mass of 2.3 nuclei per hydrogen molecule), is different from 
the sound speed in the envelope, Cs, because th e midplane tem- 
perature of the disk, T^, varies as described in Hueso & GuillotI 
(2005). The other variable from Eq. (fT4l i is the scale height: 



where Ok is the Keplerian rotation rate: 



GM, 



(15) 



(16) 



with M» the stellar mass [Eq. (|2 9]l]. The viscosi t y pa- 
rameter a is kept constant at 10'^ dPuUemond et alj uOOTI : 
[Andrews & Winiamsll2007bl) . 

Solving the problem of advectio n and diffusion yields 
the radial velocities inside the disk dPuUemond et alj I2OO6I: 
iLvnden-BeU & Pringlel[T974l) : 



URiRj)^ p— (SvV^) . 

•£.^|R^R^ ' 

The surface density evolves as 



(17) 



1 Q 

--—{I.Rur) + S 
RdR 



(18) 



dI.(R, t) 
dt 

where the source function S accounts for the infall of material 
from the envelope: 



S{R,t) = 2Npu,, 



(19) 



with u, the vertical component of the envelope velocity field 
[Eqs. dSI, dTOl ) and dTT])1. The factor 2 accounts for the enve- 
lope accreting onto both sides of the disk and the normalization 
factor ensures that the overall accretion rate onto the star and 
the disk is always equal to M. Both p and u, in Eq. dT9] l are to be 
computed at the disk-envelope boundary, which will be defined 

at the end of thi s section. 

As noted bv lHueso & Guillol d2005h . the infalling envelope 
material enters the disk with a subkeplerian rotation rate, so, 
by conservation of angular momentum, it would tend to move 
a bit further inwards. Not taking this into account would arti- 
ficially generate angular momentum, causing the disk to take 
longer to accrete onto the star As a consequence the disk will, 
at any given point in tim e, have too high a mass and too large a 
radius. iHueso & Guillot solved this problem by modifying Eq. 
( fT9] l to place the material directly at the correct radius. However, 
this causes an undesirable discontinuity in the infall trajectories. 
Instead, our model adds a small extra component to Eq. dTTl l for 



t < f.„ 



, . 3 5 / /-N [GM 



(20) 
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The functional form of the extra term derives from the CM81 
solution. A constant value of 0.002 for rji- is found to reproduce 
very well the results of Hueso & Guillot. It also provides a goo d 
match wit h the disk masses from lYorke & Bodenheimerl (1 19991) . 
lYE05l and lBWH08L whose models cover a wide range of initial 
conditions. 

The disk's inner radius is determined by the evaporation of 
dust by the star (e.g. .YE05i) : 



Ri(t) 



\ 47rcrre'tap 



(21) 



where cr is the Stefan-Boltzmann constant. The dust evaporation 
temperature, Tevap, is set to 2000 K. Taking an alternative value 
of 1500 K has no effect on our results. The stellar luminosity, 
L,, is discussed in Section l23] Inward transport of material at 
leads to accretion from the disk onto the star: 



(22) 



with the radial velocity, mr, and the surface density, 2, taken at 
Ri. The disk gains mass from the envelope at a rate Mg^d, so the 
disk mass evolves as 

Md(f) = j| (Me^d - Md^.) df' . (23) 

Our model us es a Gaussian profile for the vertical structure 
of the disk (.Shakura & Sunvae' 



(24) 



ipronle i 
vll973l) : 



p(7?,z,f) =Pcexp(-— 1 , 



with z the height above the midplane. The scale height comes 
from Eq. ( fTSl ) and the midplane density is 



(25) 



Along with the radial motion [Eq. ( l20l i. taken to be indepen- 
dent of z], material also moves vertically in the disk, as it must 
maintain the Gaussian profile at all times. To see this, consider a 
parcel of material that enters the disk at time t at coordinates R 
and z into a column with scale height H and surface density E. 
The column of material below the parcel is 



rp(/?,^,f)d^=^2:erf 
Jo 2 \//V2 



(26) 



where erf is the error function. At a later time f', the entire col- 
umn has moved to R' and has a scale height H' and a surface 
density S'. The same amount of material must still be below the 
parcel: 



Rearranging gives the new height of the parcel, z': 



z!{R!X) = //' V2erf" 



— erf 

^' \//V2 



(27) 



(28) 



where erf is the inverse of the error function. In the absence of 
mixing, our description leads to purely laminar flow. 

The location of the disk-envelope boundary [needed, for in- 
stance, in Eq. ( fT9] l1 is determined in two steps. First, the surface 
is identified where the density due to the disk [Eq. ( I24l l1 equals 




Fig. 1. Schematic view of the disk-envelope boundary in the 
(R,z) plane. The black line indicates the surface where the den- 
sity due to the disk equals that due to the envelope. The grey line 
is the infall trajectory that would lead to point Pi. However, it 
already intersects the disk at point P^, so no accretion is possi- 
ble at Pi. The disk-envelope boundary is therefore raised at Pi 
until it can be reached freely by an infall trajectory. 

that due to the envelope [Eqs. Q and (fT2li1. In order for accre- 
tion to take place at a given point Pi on the surface, it must 
be intersected by an infall trajectory. Due to the geometry of 
the surface, such a trajectory might also intersect the disk at a 
larger radius P2 (Fig. [TJ. Material flowing in along that trajec- 
tory will accrete at P2 instead of Pi. Hence, the second step in 
determining the disk-envelope boundary consists of raising the 
surface at "obstructed points" like Pi to an altitude where ac- 
cretion can take place. The source function is then computed at 
that altitude. Physically, this can be understood as follows: the 
region directly above the obstructed points becomes less dense 
than what it would be in the absence of a disk, because the disk 
also prevents material from reaching there. The lower density 
above the disk reduces the downward pressure, so the disk puffs 
up and the disk-envelope boundary moves to a higher altitude. 

The infall trajectories in the vicinity of the disk are very 
shallow, so the bulk of the material accretes at the outer edge. 
Because the disk quickly spreads beyond the centrifugal radius, 
much of the accretion occurs far from the star. In contrast, ac- 
cretion in one-dimensional collapse models occurs at or inside 
of R^ . Our r esults are consistent with the hydrodynamical work 
of BWH08, where most of the infalling material also hits the 
outer edge of a rather large disk. The large accretion radii lead to 
weaker accretion shocks than commonly assumed (Section |23]l. 

2.3. Sfar 

The star gains material from the envelope and from the disk, so 
its mass evolves as 



M,(f) = ^ (Me^. +Md^.)df'. 



(29) 



The protostar does not come into existence immediately at the 
onset of collapse; it is preceeded by the first hyd rostatic core 
(FHC; lMasunagaet"anil998l: iBoss & York3ll995h . Our model 
follows YE05 and takes a lifetime of 2 x 10^ yr and a size of 
5 AU for the FHC, independent of other parameters. After this 
stage, a rapid transition occurs from the large FHC to a protostar 
of a few P©: 



P» = (5AU) 



1 - 



f-20, 000 yr 



\ 100 yr 
20,000 <f(yr)< 20, 100, 



+ R 



PS 



(30) 
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Fig. 2. Evolution of the mass of the envelope, star and disk (left 
panel) and the luminosity (solid lines) and radius (dotted lines) 
of the star (right panel) for our standard model (black lines) and 
our reference model (grey lines). 



wher e RZ^ (ranging from 2 to 5 Rq) is the protostellar radius 
from lPalla & Stahled(ll991l) . For t > 2.01 x lO^ yr, R, equals R^^. 
Our results are not sensitive to the exact values used for the size 
and Ufetime of the FHC or the duration of the FHC-protostar 
transition. 

The star's luminosity, L,, consists of two terms: the accre- 
tion luminosity, Lt acc, dominant at early times, and the luminos- 
ity due to gravitational contraction an d deuterium burning, Z ^hot- 
The accretion luminosity comes from lAdams & Shul (Il986h : 



i.,acc(0 = '^('1^ [6m, - 2 + (2 - 5m.) Vl -m.] + 

^ [1 - (1 - 77d)Md] [l - (1 - Jld) Vl - «»] I > (31) 

where Lq - GMM/Rt. (with M the total accreted mass, i.e., M - 
M« + Md), M, = Rt/Rc, and 



Md 



du . 



(32) 



Analytical solutions exist for the asymptotic cases of m» ^0 
and M, =s 1 . For intermediate values, the integral must be solved 
numerically. The efficiency parameters 77, and rj^ in Eq. ( |3TT l 
have values of 0.5 and 0.75 for a 1 M r:, envelope (,YE05). The 
photos pheric lumino sity is adopted from lD'Antona & MazzitelUl 
using'YE05''s method of fitting and interpolating, includ- 
ing a time difference of 0.38 face (e qual to the free-fal l time) be- 
tween the onset of L, and L,,phot (iMvers et al.ll998l) . The sum 
of these two terms gives the total stellar luminosity: 



phot ' 



(33) 



Figure |2] shows the evolution of the stellar mass, luminosity 
and radius for our standard case of Mq = 1 .0 Mq, Cs = 0.26 km 
s-i and Qo = lO"''* s-\ and our reference case of Mq - 1.0 
Mo, Cs = 0.26 km s"' and Qq = 10"'^ s^' (Section |221i. The 
transition from the FHC to the protostar at t - 2 x 10'* yr is 
clearly visible in the and L, profiles. At f = face, there is no 
more accretion from the envelope onto the star, so the luminosity 
decreases sharply. 

The masses of the disk and the envelope are also shown in 
Fig. |2] Our disk mass of 0.43 Mq at f = face in the reference 
case is in excellent agreement with the value of 0.4 Mq found by 
iBWHOSl forthe same parameters. 




50 100 150 200 



20 40 60 80 100 
(K) 



Fig. 3. Dust temperature due to the accretion shock (vertical 
axis) and stellar radiation (horizontal axis) at the point of entry 
into the disk for 0.1 -//m grains in a sample of several hundred 
parcels in our standard (left) and reference (right) models. These 
parcels occupy positions from R - I to 300 AU in the disk at 
face- Shock heating is not important for the dust in most of these 
parcels, because the shock temperature is lower than the radia- 
tive heating temperature for all of them. Note the different scales 
between the two panels. 



2.4. Temperature 

The envelope starts out as an isothermal sphere at 10 K 
and it is heated up from the inside after the onset of col- 
lapse. Using the star as the only photon source, the dust tem- 
perature in the disk and envelope is computed with the ax- 
isymmetric three-dimensional radiative transfer code RADMC 
dPuUemond & DominikI l2004 l). Because of the high densities 
throughout most of the system, the gas and dust are expected to 
be well coupled, and the gas temperature is set equal to the dust 
temperature. Cosmic -ray heating of the gas is included implic- 
itly by setting a lower limit of 8 K in the dust radiative transfer 
results. As mentioned in Section lTTl the presence of the outflow 
cones has some effect on the temperature profiles ( Whitney et al.l 
I2OO3I) . This will be discussed further in Section [J!2] 



2.5. Accretion shock 

The infall of high-velocity envelope material into the low- 
velocity disk causes a J-type shock. The temperature right be- 
hind the shock front can be much higher than what it would 
be due to the stellar photons. Neufeld & Hollenbach ( 1994) cal- 
culated in detail the relationship between the pre-shock veloci- 
ties and densities (Ms and «s) and the maximum grain tempera- 
ture reached after the shock {Td,s)- A simple formula, valid for 
Ms < 70 km s can be extracted from their Fig. 13: 

'^■^^^^''Mwt^) (30^)(aT7^) '^^^^ 

with flgr the grain radius. The exponent p is 0.62 for Ms < 30 km 
s"' and 1.0 otherwise. 

The pre-shock velocities and densities are highest at early 
times, when accretion occurs close to the star and all ices would 
evaporate anyway. Important for our purposes is the question 
whether the dust temperature due to the shock exceeds that due 
to stellar heating. If all grains have a radius of 0.1 fim, as as- 
sumed in our model, this is not the case for any of the material 
in the disk at face for either our s tandard or our reference model 
(Fig. [2 cf. lSimonelli e"tanfT997h . 
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In reality, the dust spans a range of sizes, extending down to 
a radius of about 0.005 fim. Small grains are heated more eas- 
ily; 0.005-/im dust reaches a shock temperature almost twice as 
high as does 0. 1 -yum dust [Eq. ( [34l )1 . This is enough for the shock 
temperature to exceed the radiative heating temperature in part 
of the sample in Fig. [3] However, this has no effect on the CO 
and H2O gas-ice ratios. In those parcels where shock heating 
becomes important for small grains, the temperature from radia- 
tive heating lies already above the CO evaporation temperature 
of about 1 8 K and the shock temperature remains below 60 K, 
which is not enough for H2O to evaporate. Hence, shock heating 
is not included in our model. 

H2O may also be removed from the grain surfac es in the ac- 
cretio n shock through sputtering (Tielens et al. 1994: lJones et alj 
1 1994 ). The material that makes up the disk at the end of the col- 
lapse in our standard model experiences a shock of at most 8 km 
s"^ At that velocity, He^, the most important ion for sputtering, 
carries an energy of 1.3 eV. However, a minimu m of 2.2 eV is 
required to remove H2O (iBohdanskvetalJI 19801) . so sputtermg 
is unimportant for our purposes. 

Some of the material in our model is heated to more than 100 
K during the collapse (Fig. [T2b or experiences a shock strong 
enough to induce sputtering. This material normally ends up in 
the star before the end of the collapse, but mixing may keep some 
of it in the disk. The possible consequences are discussed briefly 
in Sectionl44l 



Table 1. Summary of the parameter grid used in our model.' 



2.6. Model parameters 

The standa rd set of parameters for our model corresponds to 
Case J from lYorke & Bodenheimerld 19991) . except that the solid- 
body rotation rate is reduced from 10"'"' to lO"'"' s"' to produce 
a more re alistic disk mass of 0.05 M q, consistent with observa- 
tions (e.g. [Andrews & Williamsir2007a b). The envelope has an 
initial mass of 1 .0 Mq and a radius of 6700 AU, and the effective 
sound speed is 0.26 km s"' . 

The original Case J (with Qq = 10 s '), which was also 
used in 'BWH08', is used here as a reference model to enable a 
direct quantitative comparison of the results with an indepen- 
dent method. This case results in a much higher disk mass of 
0.43 Mq. Although such high disk mass es are not excluded b y 
observations and theoretical arguments jHartmann et al.ll2006h . 
they are considered less representative of typical young stellar 
objects than the disks of lower mass. 

The parameters Mq, Cs and Qq are changed in one direc- 
tion each to create a 2^ parameter grid. The two values for Qo, 
10"'^ and 10~'^ s~', co ver the range of rotation rates observed by 
iGoodman et alj (Il993b . The other variations are chosen for their 
opposite effect: a lower sound speed gives a more massive disk, 
and a lower initial mass gives a less massive disk. The full model 
is run for each set of parameters to analyse how the chemistry 
can vary between different objects. The parameter grid is sum- 
marised in Table [1] Our standard set is Case 3 and our reference 
set is Case 7. 

Table [T] also lists the accretion time and the adsorption 
timescal e for H2O at the ed ge of the initial envelope. For com- 
parison, lEvans et alj (l2009h found a median timescale for the 
embedded phase (Class and I) of 5.4 x 10^ yr from observa- 
tions. It should be noted that the end point of our model (face) 
is not yet representative of a typical T Tauri disk (see Section 
I3.2I 1. Nevertheless, it allows an exploration of how the chemistry 
responds to plausible changes in the environment. 



Case^ 


Do 




Mo 


4cc 


Tads 


Mi 




(s ') 


/I 1 \ 

(km s ') 


(Me) 


(10^ yr) 


(10'' yr) 


(Mo) 


1 


ID 


n 1 n 


1 n 
l.U 


o.i 


1/1/1 
14.4 


U./Z 


2 


10-" 


0.19 


0.5 


3.2 


3.6 


0.08 


3 (std) 


10-'" 


0.26 


1.0 


2.5 


2.3 


0.05 


4 


10-'* 


0.26 


0.5 


1.3 


0.6 


0.001 


5 


10-'3 


0.19 


1.0 


6.3 


14.4 


0.59 


6 


10-'3 


0.19 


0.5 


3.2 


3.6 


0.25 


7 (ref) 


io-'3 


0.26 


1.0 


2.5 


2.3 


0.43 


8 


io-'3 


0.26 


0.5 


1.3 


0.6 


0.16 



'' Qq: solid-body rotation rate; c^: effective sound speed; Mq: initial 
envelope mass; t^^c- accretion time; r^ds: adsorption timescale for 
H2O at the edge of the initial cloud; M^: disk mass at face-. 
Case 3 is our standard parameter set and Case 7 is our reference set. 



2. 7. Adsorption and desorption 

The adsorption and desorption of CO and H2O are solved in a 
Lagrangian frame. When the time-dependent density, velocity 
and temperature profiles have been calculated, the envelope is 
populated by a number of parcels of material (typically 12,000) 
at f = 0. They fall in towards the star or disk according to the 
velocity profiles. The density and temperature along each par- 
cel's infall trajectory are used as input to solve the adsorption- 
desorption balance. Both species start fully in the gas phase. 
The envelope is kept static for 3 x 10^ yr before the onset of 
collapse to simulate the pre-stellar core phase. This is the same 
value as used by Rodgers & Charnlev (2003) and BWH08, and 
it is consistent with recent observations bv Enoch et al.l (l2008h . 
The amount of gaseous material left over near the end of the 
pre-collapse phase is also consistent with observations, which 
show that the onset of H2O ice formation is around an Ay of 3 
(Whittet et al. 2001). In six of our eight parameter sets, the ad- 
sorption timescales of H2O at the edge of the cloud are shorter 
than the combined collapse and pre-collapse time (Table [Til, 
so all H2O is expected to freeze out before entering the disk. 
Because of the larger cloud size, the adsorption timescales are 
much longer in Cases 1 and 5, and some H2O may still be in the 
gas phase when it reaches the disk. 

No chemical reactions are included other than adsorption and 
thermal desorption, so the total abundance of CO and H2O in 
each parcel remains constan t. Th e adsorption rate in cm-^s-' is 
taken from iCharnlev et al.i (i2001h : 



R-adsiX) = (4.55 X 10"'^ cm^ K'^'^ s-')«H«g(X) - 



To 



(35) 



where rtu is the total hydrogen density, Tg the gas temperature, 
«g(X) the gas-phase abundance of species X and M{X) its molec- 
ular weight. The numerical factor assumes unit sticking effi- 
ciency, a grain radius of 0. 1 fim and a grain abundance «gi- of 
10 '^ with respect to H2. 

The thermal desorption of CO and H2O is a zero th-order pro- 
cess: 



^des(^) = (1-26 X 10"^' cm^)nH/(^)v(X) exp 



where Td is the dust temperature and 
ns(X)] 



Ei(X) 



fiX) = min 



1, 



Nbn„ 



(36) 



(37) 
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Table 2. Binding energies and desorbing fractions for the four- 
flavour CO evaporation model. ^ 



Flavour 


£b(CO)A (K)" 


Fraction"^ 


1 


855 


0.350 


2 


960 


0.455 


3 


3260 


0.130 


4 


5773 


0.065 



^ Based on lWti et al.l ( l2004h . 
The rates for Flavours 1-3 are computed from Eq. l |36t with X = 
CO. The rate for Flavour 4 is equal to the H2O desorption rate. 
These numbers indicate fractions of adsorbing CO: 35% of all ad- 
sorbing CO becomes Flavour 1, and so on. 



with n^{X) the solid abundance of species X and A^b - 10^ the 
typical number of binding sites per grain. The numerical factor 
in Eq. ( |36] | assumes the same grain properties as in Eq. (|35] |. The 
pre-exponential factor, v{X), and the binding energy, Eb(X)/k, 
are set to 7 x 10-'' cm^^ g-i ^jjjj §55 fQj. cp and to 1 x 10^° 
cm'^ s^' and 5773 K for H2O (iBisschop et al.ll2006l:lFraser et all 
I2001h . 

Using a single Eb{CO) value means that all CO evaporates 
at the same temperature. This would be appropriate for a pure 
CO ice, but not for a mixed CO-H2O ice as is likely to form 
in reality. During the warm-up phase, part of the CO is trapped 
inside the H2O ice until the temperature becomes high enough 
for the H2O to evaporate. Recent laboratory experiments suggest 
that CO desorbs from a CO-H2O ice in four steps dCoUings et alj 
l2004i) . This can be simulated with four "flavours" of CO ice, 
each wifli a diff'erent Eb{CO) value dViti et alj|2004l) . For each 
flavour, the desorption is assumed to be zeroth order. The four- 
flavour model is summarised in Table |2l 



3. Results 

Results are presented in this section for our standard and refer- 
ence models (Cases 3 and 7) as described in Section [231 These 
cases will be compared to the other parameter sets in Section l4n 
Appendix |A] describes a formula to estimate the disk formation 
efficiency, defined as Mi/Mo at the end of the collapse phase, 
based on a fit to our model results. 

3.1. Density profiles and infall trajectories 

In our standard model (Case 3), the envelope collapses in 2.5 x 
10^ yr to give a star of 0.94 Mq and a disk of 0.05 Mq. The re- 
maining 0.01 Mq has disappeared through the bipolar outflow. 
The centrifugal radius in our standard model at face is 4.9 AU, 
but the disk has spread to 400 AU at that time due to angular mo- 
mentum redistribution. The densities in the disk are high: more 
flian 10^ cm"^ at the midplane inside of 120 AU (Fig. |4] left) 
and more than 10'^ cm"^ near 0.3 AU. The corresponding sur- 
face densities of the disk are 2.0 g cm-^ at 120 AU and 660 g 
cm-" at 0.3 AU. 

Due to the higher rotation rate, our reference model (Case 7) 
gets a much higher disk mass: 0.43 Mq. This value is consistent 
with the mass of 0.4 Mq reported by BWH08. Overall, the ref- 
erence densities from our semi-analytical model (Fig. |4] right) 
compare well with those from their more realistic hydrodynam- 
ical simulations; the dififerences are generally less than a factor 
of two. 



In both cases, the disk first emerges at 2 x 10 yr, when the 
FHC contracts to become the protostar, but it is not until a few 
lO'* yr later that the disk becomes visible on the scale of Fig. 
|4] The regions of high density (nu > 10^-10'' cm"^) are still 
contracting at that time, but the growing disks eventually cause 
them to expand again. 

Material falls in along nearly radial streamlines far out in the 
envelope and deflects towards the midplane closer in. When a 
parcel enters the disk, it follows the radial motion caused by the 
viscous evolution and accretion of more material from the en- 
velope. At any time, conservation of angular momentum causes 
part of the disk to move inwards and part of it to move outwards. 
An individual parcel entering the disk may move out for some 
time before going further in. This leads the parcel through sev- 
eral density and temperature regimes, which may aff'ect the gas- 
ice ratios or the chemistry in general. The back-and-forth motion 
occurs especially at early times, when the entire system changes 
more rapidly than at later times. The parcel motions are visu- 
alised in Figs. |5] and |6] where infall trajectories are drawn for 24 
parcels ending up at one of eight positions at facc^ at the midplane 
or near the surface at radial distances of 10, 30, 100 and 300 AU. 
Only parcels entering the disk before f ^ 2 x 10^ yr in our stan- 
dard model or f =s 1 x 10^ yr in our reference model undergo the 
back-and-forth motion. The parcels ending up near the midplane 
all entered the disk earlier than the ones ending up at the surface. 

Accretion from the envelope onto the disk occurs in an 
inside-out fashion. Because of the geometry of the disk (Fig.[TJ, 
most of the material enters near the outer edge and prevents the 
older material from moving further out. Our flow inside the disk 
is purely laminar, so some material near the midplane does move 
outwards underneath the newer material at higher altitudes. 

Because of the low rotation rate in our standard model, the 
disk does not really begin to build up until 1.5 x 10^ yr (0.6racc) 
after the onset of collapse. In addition, most of the early material 
to reach the disk proceeds onto the star before the end of the 
accretion phase, so the disk at face consists only of material from 
the edge of the original cloud (Fig.|7] left two panels). 

The disk in our reference model, however, begins to form 
right after the FHC-protostar transition at 2 x 10"^ yr. As in the 
standard model, a layered structure is visible in the disk, but it is 
more pronounced here. At the end of the collapse, the midplane 
consists mostly of material that was originally close to the centre 
of the envelope (Fig.|7] right two panels). The surface and outer 
parts of the disk are made up primarily of material from the outer 
parts of the envelope. This was also reported bv lBWHOSl 

3.2. Temperature profiles 

When the star turns on at 2 x 10"^ yr, the envelope quickly heats 
up and reaches more than 100 K inside of 10 AU. As the disk 
grows, its interior is shielded from direct irradiation and the mid- 
plane cools down again. At the same time, the remnant envelope 
mat erial above the di sk becomes less dense and warmer As in 
Whitney et al.l (l2003h . the outflow has some effect on the tem- 
perature profile. Photons emitted into the outflow can scatter and 
illuminate the disk from the top, causing a higher disk tempera- 
ture beyond R x 200 AU than if there were no outflow cone. At 
smaller radii, the disk temperature is lower than in a no-outflow 
model. Without the outflow, the radiation would be trapped in 
the inner envelope and inner disk, increasing the temperature at 
small radii. 

At f = face in our standard model, the 100- and 18-K 
isotherms (where H2O and pure CO evaporate) intersect the mid- 
plane at 20 and 2000 AU (Fig. [8] left). The disk in our reference 
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Fig. 4. Total density at four time steps for our standard model (Case 3; left) and our reference model (Case 7; right). The time is 
given in years as well as in units of the accretion time, with a{b) meaning a x 10^. The density contours increase by factors of ten 
going inwards; the 10^-cm"^ contours are labelled in the standard panels and the 10''-cm"^ contours in the reference panels. The 
white curves indicate the surface of the disk as defined in Section l2!2l (onlv visible in three panels). Note the different scale between 
the two sets of panels. 



model is denser and therefore colder: it reaches 100 and 18 K 
at 5 and 580 AU on the midplane (Fig. [8] right). Our radiative 
transfer method is a more rigorous way to obtain t he dust te m- 
perature than the diffusion approximation used by lBWH08l so 
our temperature profiles are more realistic than theirs. 

Compared to typ ical T Tau ri disk models (e.g. 

iD'Alessio etakl 119981 1 19991 1200 ll) . our standard disk at 
face is warmer. It is 81 K at 30 AU on the midplane, while the 
closest model from the D'Alessio catalogue is 28 K at that 
point. If our model is allowed to run beyond face, part of the disk 
accretes further onto the star. At f = 4facc (10^ yr), the disk mass 
goes down to 0.03 Mq. The luminosity of the s tar decreases 
during this period (ID Antona & Mazzitellil Il994t) . so the disk 
cools down: the midplane temperature at 30 AU is now 42 K. 
Meanwhile, the dust is likely to grow to larger sizes, whic h 
would further decrease the temperatures (I D'Alessio et al.ll2001h . 
Hence, it is important to realise that the normal end point of our 
models does not represent a "mature" T Tauri star and disk as 
typically discussed in the literature. 



3.3. Gas and ice abundances 

Our two species, CO and H2O, begin entirely in the gas phase. 
They freeze out during the static pre-stellar core phase from the 
centre outwards due to the density dependence of Eq. (l35b . After 
the pre-collapse phase of 3 x 10^ yr, only a few tenths of per cent 
of each species is still in the gas phase at 3000 AU. About 30% 
remains in the gas phase at the edge of the envelope. 

Up to the point where the disk becomes important and the 
system loses its spherical symmetry, our model gives the same 
results as the one-dimensional collapse models: the temperature 
quickly rises to a few tens of K in the collapsing region, driving 
some CO (evaporating around 18 K in the one-flavour model) 
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Fig. 9. Gaseous CO as a fraction of the total CO abundance (top) 
and idem for H2O (bottom) at two time steps for our standard 
model (Case 3). The black curves indicate the surface of the disk 
(only visible in two panels). The black area near the pole is the 
outflow, where no abundances are computed. Note the different 
spatial scale between the two panels of each set; the small box 
in the left CO panel indicates the scale of the H2O panels. 
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Fig. 5. Infall trajectories for parcels in our standard model (Case 3) ending up near the surface (top panels) or at the midplane 
(bottom panels) at radial positions of 10, 30, 100 and 300 AU (dotted lines) at f = t.^^^. Each panel contains trajectories for three 
parcels, ■which are illustrative for material ending up at the given location. Trajectories are only dra'wn up to f = face- Diamonds 
indicate ■where each parcel enters the disk; the time of entry is given in units of 10^ yr. Note the different scales between some 
panels. 
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Fig. 6. Same as Fig. |5] but for our reference model (Case 7). 



into the gas phase, but keeping H2O (evaporating around 100 K) 
on the grains. 

As the disk grows in mass, it provides an increasingly 
large body of material that is shielded from the star's radia- 
tion, and that is thus much colder than the surrounding envelope. 
However, the disk in our standard model never gets below 1 8 K 
before the end of the collapse (Section lTSb . so CO remains in the 
gas phase (Fig.|9] top). Note that trapping of CO in the H2O ice 



is not taken into account here; this possibility will be discussed 
in Sectionl43] 



The disk in our reference model is more massive and there- 
fore colder After about 5x10^ yr, the outer part drops below 18 
K. CO arriving in this region re-adsorbs onto the grains (Fig.fTOl 
top). Another 2 x 10^ yr later, at f = f^cc, 19% of all CO in the 
disk is in solid form. Moving out from the star, the first CO ice 
is found at the midplane at 400 AU. The solid fraction gradually 
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Standard: Case 3 Reference: Case 7 




Fig. 7. Position of parcels of material in our standard model (Case 3; left) and our reference model (Case 7; right) at the onset of 
collapse (t - 0) and at the end of the collapse phase (f - t^^c)- The parcels are colour-coded according to their initial position. 
A layered structure is visible in the disk, with material near the surface and at the outer edge originating from further out in the 
envelope than the material near the midplane. This effect is most pronounced in our reference model. The grey parcels from t - 
are in the star or have disappeared through the outflow at f = t^^c- Note the different spatial scale between the two panels of each set; 
the small box in the left panel indicates the scale of the right panel. 
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Fig. 8. Dust temperature, as in Fig.El Contours are drawn at 100, 60, 50, 40, 35, 30, 25, 20, 18, 16, 14 and 12 K. The 40- and 20-K 
contours are labelled in the standard and reference panels, respectively. The 18- and 100-K contours are drawn as thick grey lines. 
The white curves indicate the surface of the disk as defined in Section l272l (onlv visible in four panels). 



increases to unity at 600 AU. At R - 1000 AU, nearly aU CO 
is solid up to an altitude of 170 AU. The solid and gaseous CO 
regions meet close to the 18-K surface. The densities throughout 
most of the disk are high enough that once a parcel of material 
goes below the CO desorption temperature, all CO rapidly dis- 
appears from the gas. The exception to this rule occurs at the 
outer edge, near 1500 AU, where the adsorption and desorption 
timescales are longer than the dynamical timescales of the in- 
falling material. Small differences between the trajectories of 
individual parcels then cause some irregularities in the gas-ice 
profile. 

The region containing gaseous H2O is small at all times dur- 
ing the collapse. At f = face, the snow line (the transition of H2O 
from gas to ice) lies at 15 AU at the midplane in our standard 



model (Fig. |9l bottom). The surface of the disk holds gaseous 
H2O out to 7? = 41 AU, and overafl 13% of afl H2O in the disk is 
in the gas phase. This number is much lower in the colder disk 
of our reference model: only 0.4%. The snow line now lies at 7 
AU and gaseous H2O can be found out to 17 AU in the disk's 
surface layers (Fig.[l0l bottom). 

Using the adsorption-desorption history of all the individual 
infalling parcels, the original envelope can be divided into sev- 
eral chemical zones. This is trivial for our standard model. All 
CO in the disk is in the gas phase and it has the same qualitative 
history: it freezes out before the onset of collapse and quickly 
evaporates as it falls in. H2O also freezes out initially and only 
returns to the gas phase if it reaches the inner disk. 
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Reference: Case 7: CO 




Reference: Case 7: H^O 
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Fig. 10. Same as Fig. |9] but for our reference model (Case 7). 
The CO gas fraction is plotted on a larger scale and at two addi- 
tional time steps. 



Our reference model has the same general H2O adsorption- 
desorption history, but it shows more variation for CO, as illus- 
trated in Fig.[TT] For the red parcels in that figure, more than half 
of the CO always remains on the grains after the initial freeze- 
out phase. On the other hand, more than half of the CO comes 
off the grains during the collapse for the green parcels, but it 
freezes out again inside the disk. The pink parcels, ending up 
in the inner disk or in the upper layers, remain warm enough to 
keep CO off the grains once it first evaporates. The blue parcels 
follow a more erratic temperature profile, with CO evaporating, 
re-adsorbing and evaporating a second time. This is related to 
the back-and-forth motion of some material in the disk (Fig.|6]l. 

3.4. Temperature histories 

The proximity of the CO and H2O gas-ice boundaries to the 18- 
and 100-K surfaces indicates that the temperature is primarily 
responsible for the adsorption and desorption. At nn - 10'' cm"^, 
adsorption and desorption of CO are equally fast at Td = 18 K 
(a timescale of 9 X 10^ yr). For a density a thousand times higher 
or lower, the dust temperature only has to increase or decrease 
by 2-3 K to maintain k^d^ = ^des- 

The exponential temperature dependence in the desorption 
rate [Eq. ( [36] l1 also holds for other species than CO and H2O, 
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Fig. 11. Same as Fig.|7] but only for our reference model (Case 7) 
and with a different colour scheme to denote the CO adsorption- 
desorption behaviour In all parcels, CO adsorbs during the 
pre-collapse phase. Red parcels: CO remains adsorbed; green 
parcels: CO desorbs and re-adsorbs; pink parcels: CO desorbs 
and remains desorbed; blue parcels: CO desorbs, re-adsorbs and 
desorbs once more. The fraction of gaseous CO in each type of 
parcel as a function of time is indicated schematically in the in- 
set in the right panel. The grey parcels from f = are in the star 
or have disappeared through the outflow at f = face- In our stan- 
dard model (Case 3), all CO in the disk at face is in the gas phase 
and it all has the same qualitative adsorption-desorption history, 
equivalent to the pink parcels. 

as well as for the rates of some chemical reactions. Hence, it is 
useful to compute the temperature history for infalling parcels 
that occupy a certain position at face. Figures [T2] and IATT3] (the 
latter is only available online) show these histories for material 
ending up at the midplane or near the surface of the disk at radial 
distances of 10, 30, 100 and 300 AU. Parcels ending up inside of 
10 AU have a very similar temperature history as those ending 
up at 10 AU, except that the final temperature of the former is 
higher. 

Each panel in Figs. [T2l and IA.13l contains the history of sev- 
eral dozen parcels ending up close to the desired position. The 
qualitative features are the same for all parcels. The temperature 
is low while a parcel remains far out in the envelope. As it falls 
in with an ever larger velocity, there is a temperature spike as 
it traverses the inner envelope, followed by a quick drop once 
it enters the disk. Inward radial motion then leads to a second 
temperature rise; because of the proximity to the star, this one is 
higher than the first increase. For most parcels in Figs. [T2l and 
IA.13I the second temperature peak does not occur until long af- 
ter face- In all cases, the shock encountered upon entering the 
disk is weak enough that it does not heat the dust to above the 
temperature caused by the stellar photons (Fig.[3]l. 

Based on the temperature histories, the gas-ice transition at 
the midplane would lie inside of 10 AU for H2O and beyond 300 
AU for CO in both our models. This is indeed where they were 
found to be in Section 13.31 The transition for a species with an 
intermediate binding energy, such as H2CO, is then expected to 
be between 10 and 100 AU, if its abundance can be assumed 
constant throughout the collapse. 

The dynamical timescales for the infalling material before it 
enters the disk are between 10^ and 10^ yr. The timescales de- 
crease as it approaches the disk, due to the rapidly increasing ve- 
locities. Once inside the disk, the material slows down again and 
the dynamical timescales return to 10'*-10^ yr The adsorption 
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Fig. 12. Temperature history for parcels in our standard model (Case 3) ending up near the surface (top panels) or at the midplane 
(bottom panels) at radial positions of 10, 30, 100 and 300 AU at f = face- Each panel contains between 22 and 90 curves; the coloured 
curves correspond to the parcels from Fig.|5] The first peak in each curve occurs just before that parcel enters the disk. The dotted 
lines are drawn at Td = 18 K and f = face- Note the different vertical scales between some panels. 



timescales for CO and H2O are initially a few 10^ yr, so they ex- 
ceed the dynamical timescale before entering the disk. Depletion 
occurs nonetheless because of the pre-collapse phase with a du- 
ration of 3 X 10^ yr. The higher densities in the disk cause the 
adsorption timescales to drop to 100 yr or less. If the tempera- 
ture approaches (or crosses) the desorption temperature for CO 
or H2O, the corresponding desorption timescale becomes even 
shorter than the adsorption timescale. Overall, the timescales for 
these specific chemical processes (adsorption and desorption) in 
the disk are shorter by a factor of 1000 or more than the dynam- 
ical timescales. 

At some final positions, there is a wide spread in the time 
that the parcels spend at a given temperature. This is especially 
true for parcels ending up near the midplane inside of 100 AU in 
our reference model. All of the midplane parcels ending up near 
10 AU exceed 18 K during the collapse; the first one does so at 
3.5 X 10"* yr after the onset of collapse, the last one at 1.6 x 10^ 
yr. Hence, some parcels at this final position spend more than 
twice as long above 18 K than others. This does not appear to be 
relevant for the gas-ice rati o, but it is important for the formation 
of more complex species dGarrod & Herbstll2006h . This will be 
discussed in more detail in Sectionl4|2] 



4. Discussion 

4.1. Model parameters 

When the initial conditions of our model are modified (Section 
I2.6I 1, the qualitative chemistry results do not change. In Cases 3, 
4 and 8, the entire disk at face is warmer than 1 8 K, and it contains 
no solid CO. In the other cases, the disk provides a reservoir of 
relatively cold material where CO, which evaporates early on in 
the collapse, can return to the grains. H2O can only desorb in the 
inner few AU of the disk and remnant envelope. 



Figures [T4l and [TS] show the density and dust temperature at 
face for each parameter set; our standard and reference models 
are the top and bottom panel in the second column (Case 3 and 
7). Several trends are visible: 



With a lower sound speed (Cases 1, 2, 5 and 6), the overall 
accretion rate (M) is smaller so the accretion time increases 



(fa, 



). The disk can now grow larger and more massive. 



In our standard model, the disk is 0.05 Mq at face and extends 
to about 400 AU radially. Decreasing the sound speed to 0. 19 
km s ' (Case 1) results in a disk of 0.22 Mq and nearly 2000 
AU. The lower accretion rate also reduces the stellar lumi- 
nosity. These effects combine to make the disk colder in the 
low-Cs cases. 

With a lower rotation rate (Cases 1^), the infall occurs in 
a more spherically symmetric fashion. Less material is cap- 
tured in the disk, which remains smaller and less massive. 
From our reference to our standard model, the disk mass 
goes from 0.43 to 0.05 Mq and the radius from 1400 to 400 
AU. The larger accretion onto the star causes a higher lumi- 
nosity. Altogether, this makes for a small, relatively warm 
disk in the low-Qo cases. 

With a lower initial mass (Cases 2, 4, 6 and 8), there is less 
material to end up on the disk. The density profile is inde- 
pendent of the mass in a Shu-type collapse [Eq. ([U], so the 
initial mass is lowered by taking a smaller envelope radius. 
The material from the outer parts of the envelope is the last 
to accrete and is more likely, therefore, to end up in the disk. 
If the initial mass is halved relative to our standard model (as 
in Case 4), the resulting disk is only 0.001 Mq and 1 AU. Our 
reference disk goes from 0.43 Mq and 1400 AU to 0.16 Mq 
and 600 AU (Cases 7 and 8). The luminosity at face is lower 
in the high-Mo cases and the cold part of the disk (T^ < 18 
K) has a somewhat larger relative size. 
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Fig. 14. Total density at f = face for each parameter set in our grid. The numbers at the top of each panel are the parameter set number, 
the rotation rate (log Qq in s"' ), the sound speed (km s"' ) and the initial mass (Mq). The density contours increase by factors of ten 
going inwards; the 10^-cm"^ contour is labelled in each panel. The white curves indicate the surfaces of the disks; the disk for Case 
4 is too small to be visible. 
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Fig. 15. Dust temperature, as in Fig. [141 The temperature contours are drawn at 100, 60, 40, 30, 25, 20, 18, 16, 14 and 12 K from 
the centre outwards; the 20-K contour is labelled in each panel. 



iDullemond et al.l (l2006l) noted that accretion occurs closer 
to the star for a slowly rotating cloud than for a fast rotating 
cloud, resulting in a larger fraction of crystalline dust in the for- 
mer case. The same effect is seen here, but overall the accretio n 
takes place further from the star than in lDuUemond et all (l2006h . 
This is due to our taking into account the vertical structure of the 
disk. Our gaseous fractions in the Iow-Qq disks are higher than 
in the high-Qo disks (consistent with a higher crystalline frac- 



tion), but not because material enters the disk closer to the star. 
Rather, as mentioned above, the larger gas content comes from 
the higher temperatures throughout the disk. 

Combining the density and the temperature, the fractions of 
cold (T^d < 18 K), warm (7^ > 18 K) and hot (T^d > 100 K) ma- 
terial in the disk can be computed. The warm and hot fractions 
are listed in Table [3] along with the fractions of gaseous CO and 
H2O in the disk at t^^c- Across the parameter grid, 23-100% of 
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Table 3. Summary of properties at t - face for our parameter 
grid.'' 



Case 


MJM" 


J warm 


/hot'^ 


/..(CO)" 




1 


0.22 


0.69 


0.004 


0.62 


0.028 


2 


0.15 


0.94 


0.035 


0.93 


0.020 


3 (std) 


0.05 


1.00 


0.17 


1.00 


0.13 


4 


0.003 


1.00 


1.00 


1.00 


1.00 


5 


0.59 


0.15 


0.0001 


0.23 


0.11 


6 


0.50 


0.34 


0.0004 


0.27 


0.02 


7 (ref) 


0.43 


0.83 


0.003 


0.81 


0.004 


8 


0.33 


1.00 


0.028 


1.00 


0.003 



^ These results are for the one-tlavour CO desorption model. 
The fraction of the disk mass with respect to the total accreted mass 
(M = M. + Mi). 

The fractions of warm (r^ > 18 K) and hot (T^ > 100 K) mate- 
rial with respect to the entire disk. The warm fraction also includes 
material above 100 K. 

The fractions of gaseous CO and H2O with respect to the total 
amounts of CO and H2O in the disk. 



the CO is in the gas, along with 0.3-100% of the H2O. This in- 
cludes Case 4, which only has a disk of 0.0014 M©. If that one is 
omitted, at most 13% of the H2O in the disk at face is in the gas. 
The gaseous H2O fractions for Cases 1, 2, 6, 7 and 8 (at most 
a few per cent) are quite uncertain, because the model does not 
have sufficient resolution in the inner disk to resolve these small 
amounts. These fractions may be lower by up to a factor of 10 or 
higher by up to a factor of 3. 

There is good agreement between the fractions of warm ma- 
terial and gaseous CO. In Case 5, about a third of the CO gas 
at face is gas left over from the initial conditions, due to the 
long adsorption timescale for the outer part of the cloud. This 
is also the case for the majority of the gaseous H2O in Cases 
1, 5 and 6. For the other parameter sets, /hot and /gas(H20) are 
the same within the error margins. Overall, the results from the 
parameter grid show once again that the adsorption-desorption 
balance is primarily determined by the temperature, and that the 
adsorption-desorption timescales are usually shorter than the dy- 
namical timescales. 

By comparing the fraction of gaseous material at the end of 
the collapse to the fraction of material above the desorption tem- 
perature, the history of the material is disregarded. For example, 
some of the cold material was heated above 18 K during the 
collapse, and CO desorbed before re-adsorbing inside the disk. 
This may affect the CO abundance if the model is expanded to 
include a full chemical network. In that case, the results from 
Table [3] only remain valid if the CO abundance is mostly con- 
stant throughout the collapse. The same caveat holds for H2O. 

4.2. Complex organic molecules 

A full chemical network is required to analyse the gas and ice 
abundances of more complex species. While this will be a topic 
for a future paper, the current CO and H2O results, combined 
with recent other work, can already provide some insight. 

In general, the formation of organic species can be divided 
into two categories: first-generation species that are formed on 
and reside in the grain surfaces, and second-generation species 
that are formed in the warm gas phase when the first-generation 
species have evaporated. Small first-generation species (like 
CH3OH) are efficiently formed during the pre-collapse phase 



dGarrod & Herbstl|2006l) . Their gas-ice ratios should be similar 
to that of H2O, due to the similar binding energies. 

Larger first-generation species such as methyl formate 
(HCOOCH3) can be formed on the grains if material spends 
at least several lO'* yr at 20-40 K. The radicals involved in the 
surface formation of HCOOCH3 (HCO and CH3O) are not mo- 
bile enough at lower temperatures and are not formed efficiently 
enough at higher temperatures. A low surface abundance of CO 
(at temperatures above 18 K) does not hinder the formation of 
HCOOCH3: HCO and CH3O are formed from reactions of OH 
and H with H2CO, which is already formed at an earlier stage 
and which does not evaporate until ~40 K ( Garrod & Herbsll 
2006). Cosmic-ray-induced photons are available to form OH 
from H2O even i n the densest parts of the disk and envelope 
(IShen et al.ll2004l) . 

As shown in Section IT?! many of the parcels ending up near 
the midplane inside of ~300 AU in our standard model spend 
sufficient time in the required temperature regime to allow for 
efficient formation of HCOOCH3 and other complex organics. 
Once formed, these species are likely to assume the same gas- 
ice ratios as H2O and the smaller organics. They evaporate in the 
inner 10-20 AU, so in the absence of mixing, complex organics 
would only be observable in the gas phase close to the star The 
Atacama Large Millimeter/submillimeter Array (ALMA), cur- 
rently under construction, will be able to test this hypothesis. 

The gas-phase route towards complex organics involves the 
hot inner envelope (Ti > 100 K), also calle d the hot core or 
hot corino in the case of lo w-mass protostars ( ICeccareUill2004t 
iBottinelU et al.ll2004l2007h . Most of the ice evaporates here and 
a rich chemistry can take p lace if ma terial spends at least sev- 
eral 10^ yr in the hot core ( Charnlev e t al. 1992). However, the 
material in the hot inner envelope in our model is essentially in 
freefall towards the star or the inner disk, and its transit time 
of a few 100 yr is too short for compl ex organics to be formed 
abundantly (see also lSchoier et alj|200 2). Additionally, the total 
mass in this region is very small: about a per cent of the disk 
mass. In order to explain the observations of second-generation 
complex molecules, there has to be some mechanism to keep 
the material in the hot core for a longer time. Alternatively, it 
has recently been suggested that molecules typically associated 
with hot cores may in fact form on the grain surfaces as well 
dGarrod et alj|2008[) . 

4.3. Mixed CO-H2O ices 

In the results presented in Section [3] all CO was taken to des- 
orb at a single temperature. In a more realistic approach, some 
of it would be trapped in the H2O ice and desorb at higher tem- 
peratures. This was simulated with four "flavours" of CO ice, as 
summarised in Table |2] With our four-flavour model, the global 
gas-ice profiles are largely unchanged. All CO is frozen out in 
the sub- 18 K regions and it fully desorbs when the temperature 
goes above 100 K. Some 10 to 20% remains in the solid phase 
in areas of intermediate temperature. In our standard model, the 
four-flavour variety has 15% of all CO in the disk at face on the 
grains, compared to 0% in the one-flavour variety. In our refer- 
ence model, the solid fraction increases from 19 to 33%. 

The grain-surface formation of H2CO, CH3OH, HCOOCH3 
and other organics should not be very sensitive to these vari- 
ations. H2CO and CH3OH are already formed abundantly be- 
fore the onset of collapse, when the one- and four-flavour mod- 
els predict equal amounts of solid CO. H2CO is then available 
to form HCOOCH3 (via the intermediates HCO and CH3O) 
during the collapse. The higher abundance of solid CO at 20- 
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40 K in the four-flavour model could slow down the forma- 
tion of HCOOCH3 somewhat, because CO destroys the OH 
needed to form HCO (Garrod & Herbst 2006). H2CO evaporates 
around 40 K, so HCOOCH3 cannot be formed efficiently any- 
more above that temperature. On the other hand, if a multiple- 
flavour approach is also employed for H2CO, some of it remains 
solid above 40 K, and HCOOCH3 can continue to be produced. 
Overall, then, the multiple-flavour desorption model is not ex- 
pected to cause large variations in the abundances of these or- 
ganic species compared to the one-flavour model. 

4.4. Implications for comets 

Comets in our solar system are known to be abundant in CO and 
they are believed to have formed between 5 and 30 A U in the 
circumsolar disk (Bockelee-Morvanet al. 2004; Kobavas hi et al.l 
l2007h . However, the dust temperature in this region at the end of 
the collapse is much higher than 18 K for all of our parameter 
sets. This raises the question of how solid CO can be present in 
the comet-forming zone. 

One possible answer lies in the fact that even at f = face, our 
objects are still very young. As noted in Section [X21 the disks 
will cool down as they continue to evolve towards "mature" T 
Tauri systems. Given the right set of initial conditions, this may 
bring the temperature below 18 K inside of 30 AU. However, 
there are many T Tauri disk models in the literature where the 
temperature at those radii remains well above the CO evapora- 
tion temperature (e.g.'D'Al essio et al.|[l 99811200 Ih . Specifically, 
models of the minimum-mass solar nebula (MMSN) predict a 
dust temperature of -40 K at 30 AU fcecar et al.] |2006V 

A more plausible solution is to turn to mixed ices. At the tem- 
peratures computed for the comet-forming zone of the MMSN, 
10-20% of all CO may be trapped in the H2O ice. Assuming 
typical CO-H2O abundance ratios, this is entirely consistent with 
observed cometary abundances (Bockelee -Morvan et al.ll2004l) . 

Large abundance variations are possible for more complex 
species, due to the different densities and temperatures at various 
points in the comet-forming zone in our model, as well as the 
different density and temperature histories for material ending up 
at those points. This seems to be at least part of the explanation 
for the chemical diversity observed in comets. Our current model 
will be extended in a forthcoming paper to include a full gas- 
phase chemical network to analyse these variations and compare 
them against cometary abundances. 

The desorption and re-adsorption of H2O in the disk- 
envelope boundary shock has been suggested as a m ethod to trap 
noble gases in the ice an d include them in comets (lOwen et"al1 
ll992l:IOwen& Bar-Nun 1993). As shown in Sections |23] and 
13.41 a number of parcels in our standard model are heated to 
more than 100 K just prior to entering the disk. However, these 
parcels end up in the disk's surface. Material that ends up at the 
midplane, in the comet-forming zone, never gets heated above 
50 K. Vertical mixing, which is ignored in our model, may 
be able to bring the noble-gas-containing grains down into the 
comet-forming zone. 

Another option is episodic accretion, resulting in tempo- 
rary heating of the disk (Section l43] l. In the subsequent cooling 
phase, noble gases may be trapped as the ices reform. The al- 
ternative of trapping the noble gases already in the pre-collapse 
phase is unlikely. This requires all the H2O to start in the gas 
phase and then freeze out rapidly. However, in reality (contrary 
to what is assumed in our model) it is probably formed on the 
grain surfaces by hydrogenation of atomic oxygen, which would 
not allow for trapping of noble gases. 



4.5. Limitations of the model 

The physical part of our model is known to be incomplete and 
this may affect the chemical results. For exam ple, our model 
does not include radial and vertical mixing. ISemenov et al.l 
(120061) and lAikawal (l2007l) recently showed that mixing can en- 
hance the gas-phase CO abundance in the sub- 18 K regions of 
the disk. Similarly, there could be more H2O gas if mixing is 
included. This would increase the fractions of CO and H2O 
gas listed in Table [3] The gas-phase abundances can also be 
enhanced by allowing for photodesorption of the ices in addi- 
tion to the therma l desorption considered here (I Shen et aLll2004l; 
lOberg et"ani2007l 12009). Mixing and photodesorption can each 
increase the total amount of gaseous material by up to a fac- 
tor of 2. The higher gas-phase fractions are mostly found in the 
regions where the temperature is a few degrees below the des- 
orption temperature of CO or H2O. 

Accretion from the envelope onto the star and disk occurs in 
our model at a constant rate M until all of the envelope mass 
is gone. However, the lack of widespread red-shifted absorp- 
tion seen in interferometric observa tions suggests that the infall 
may stop already at an earlier time (|j0rgensen et al.ll2007l) . This 
would reduce the disk mass at face- The size of the disk is de- 
termined by the viscous evolution, which would probably not 
change much. Hence, if accretion stops or slows down before 
face, the disk would be less dense and therefore warmer It would 
also reduce the fraction of disk material where CO never des- 
orbed, because most of that material comes from the outer edge 
of the original cloud (Fig. fTTT i. Both effects would increase the 
gas-ice ratios of CO and H2O. 

Our results are a lso modified by the likely occurence 
of episodic accretion (' Ken von & HartmannI 119951 lEvans et aTl 
2009). In this scenario, material accretes from the disk onto the 
star in short bursts, separated by intervals where the disk-to-star 
accretion rate is a few orders of magnitude lower The accretion 
bursts cause luminosity flares, briefly heating up the disk before 
returning to an equilibrium temperature that is lower than in our 
models. This may produce a disk with a fairly large ice content 
for most of the time, which evaporates and re-adsorbs after each 
accretion episode. The consequences for complex organics and 
the inclusion of various species in comets are unclear. 

5. Conclusions 

This paper presents the first results from a two-dimensional, 
semi-analytical model that simulates the collapse of a molec- 
ular cloud to form a low-mass protostar and its surrounding 
disk. The model follows individual parcels of material from the 
cloud into the star or disk and also tracks their motion inside 
the disk. It computes the density and temperature at each point 
along these trajectories. The density and temperature profiles are 
used as input for a chemical code to calculate the gas and ice 
abundances for carbon monoxide (CO) and water (H2O) in each 
parcel, which are then transformed into global gas-ice profiles. 
Material ending up at different points in the disk spends a diff'er- 
ent amount of time at certain temperatures. These temperature 
histories provide a first look at the chemistry of more complex 
species. The main results from this paper are as follows: 

• Both CO and H2O freeze out towards the centre of the cloud 
before the onset of collapse. As soon as the protostar turns 
on, a fraction of the CO rapidly evaporates, while H2O re- 
mains on the grains. CO returns to the solid phase when it 
cools below 18 K inside the disk. Depending on the initial 
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conditions, this may be in a small or a large fraction of the 
disk (Section[33]). 

• All parcels that end up in the disk have the same qualitative 
temperature history (Fig.[T2b. There is one temperature peak 
just before entering the disk, when material traverses the in- 
ner envelope, and a second one (higher than the first) when 
inward radial motion brings the parcel closer to the star. In 
some cases, this results in multiple desorption and adsorption 
events during the parcel's infall history (Section [34b . 

• Material that originates near the midplane of the initial enve- 
lope remains at lower temperatures than material originating 
from closer to the poles. As a result, the chemical content 
of the material from near the midplane is less strongly mod- 
ified during the collapse than the content of material from 
other regions (Fig. [TTT i. The outer part of the disk contains 
the chemically most pristine material, where at most only a 
small fraction of the CO ever desorbed (Section [331 l. 

• A higher sound speed results in a smaller and warmer disk, 
with larger fractions of gaseous CO and H2O at the end of the 
envelope accretion. A lower rotation rate has the same effect. 
A higher initial mass results in a larger and colder disk, and 
smaller gaseous CO and H2O fractions (Section |4Tt . 

• The infalling material generally spends enough time in a 
warm zone that first-generation complex organic species can 
be formed abundantly on the grains (Fig. [T2b . Large differ- 
ences can occur in the density and temperature histories for 
material ending up at various points in the disk. These differ- 
ences allow for spatial abundance variations in the complex 
organics across the entire disk. This appears to be at least 
part of the explanation for the cometary chemical diversity 
(Sections |4!2l and |4!4l) . 

• Complex second-generation species are not formed abun- 
dantly in the warm inner envelope (the hot core or hot corino) 
in our model, due to the combined effects of the dynamical 
timescales and low mass fraction in that region (Section|4|2]). 

• The temperature in the disk's comet-forming zone (5-30 AU 
from the star) lies well above the CO desorption temperature, 
even if effects of grain growth and continued disk evolution 
are taken into account. Observed cometary CO abundances 
can be explained by mixed ices: at temperatures of several 
tens of K, as predicted for the comet-forming zone, CO can 
be trapped in the H2O ice at a relative abundance of a few 
per cent (Section |4~4] i. 
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Appendix A: Disk formation efficiency 

The results from our parameter grid can be used to derive the 
disk formation efficiency, ri^f, as a function of the sound speed, 
Cs, the solid-body rotation rate, Qo, and the initial cloud mass. 
Mo. This efficiency can be defined as the fraction of Mq that is 
in the disk at the end of the collapse phase (t = face) or as the 
mass ratio between the disk and the star at that time. The former 
is used in this Appendix. 

In order to cover a larger range of initial conditions, the phys- 
ical part of our model was run on a 9^ grid. The sound speed was 
varied from 0.15 to 0.35 km s"', the rotation rate from lO"'"*'^ to 



10 s and the initial cloud mass from 0.1 to 2.1 Mq. The 
resulting ?7df at f = face were fitted to 



log(Qo/s-i) 
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Equation lA.il can give values lower than or larger than 1. In 
those cases, it should be interpreted as being or 1 . 

The best-fit values for the coefficients kj and the exponents 
qi are listed in Table lA.ll The absolute and relative difference 
between the best fit and the model data have a root mean square 
(rms) of 0.04 and 5%. The largest absolute and relative differ- 
ence are 0.20 and 27%. The fit is worst for a high cloud mass, 
a low sound speed and an intermediate rotation rate, as well as 
for a low cloud mass, an intermediate to high sound speed and a 
high rotation rate. 

Fi gure IaTI shows the disk formation efficiency as a function 
of the rotation rate, including the fit from Eq. dA.ll ). The effi- 
ciency is roughly a quadratic function in log Q.q, but due to the 
small dynamic range of this variable, the fit appears as straight 
lines. Furthermore, the efficiency is roughly a linear function in 
Cs and a square root function in Mq. 



Table A.l. Coefficients and exponents for the best fit for the disk 
formation efficiency. 
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Fig. A.l. Disk formation efficiency as a function of the solid- 
body rotation rate. The model values are plotted as symbols and 
the fit from Eq. ( IA.lt as lines. The different values of the sound 
speed are indicated by colours and the different values of the 
initial cloud mass are indicated by symbols and line types, with 
the solid lines corresponding to the asterisks, the dotted lines to 
the diamonds and the dashed lines to the triangles. 
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Fig. A.13. Same as Fig. [12] but for our reference model (Case 7). The coloured curves correspond to the parcels from Fig. |6l 



